%%Solve for the BGP equilibrium of prices, wages, ouput, trade shares 

tic;

diff=10; iter=0; 



%...GIT'R DUN!...
% display('---------------------------------------------------------------');
 
while diff>10^-5;
     
iter=iter+1;

w=w./w(M);

Q = (TFP./Tnew).^(1/(sig-1));

TNT=(lamNT).*Tnew; %%Non-technology intensive productivity proportional to IP intensive one to guaranteee BGP where all countries and sectors grow at the same rate.

%% Compute prices

for i=1:M
    for j=1:M
PtempT(i,j)=(Tnew(j)*Q(j)^(sig-1))*(mbar*w(j)*distT(i,j)*(1+tau(i,j))).^(1-sig);
PtempNT(i,j)=(TNT(j)*Q(j)^(sig-1))*(mbar*w(j)*distNT(i,j)).^(1-sig);
    end
end
PT =sum(PtempT').^(1/(1-sig)); %%Price of T sector
PNT =sum(PtempNT').^(1/(1-sig)); %%Price NT sector

P=(PT./alpha).^alpha.*(PNT./(1-alpha)).^(1-alpha); %%Price index


%% Compute trade shares

for i=1:M
    for j=1:M
        pi_shareT(i,j)=Q(j)^(sig-1)*Tnew(j)*(distT(i,j)*(1+tau(i,j)))^(1-sig)*(w(j)*mbar)^(1-sig)/PT(i)^(1-sig); %%Trade share T sector
        pi_shareT_pretax(i,j)=pi_shareT(i,j)/(1+tau(i,j));
        pi_shareNT(i,j)=Q(j)^(sig-1)*TNT(j)*distNT(i,j)^(1-sig)*(w(j)*mbar)^(1-sig)/PNT(i)^(1-sig); %%Trade share NT sector
    end
end

%%Compute expenditures

YT=((pi_shareT_pretax)')^(-1)*(mbar.*w.*LT)./P'; %%Labor market clearing T sector
YNT=((pi_shareNT)')^(-1)*(mbar.*w.*LNT)./P'; %%Labor market clearing NT sector

Y=YT; %%Final production
Yw=sum(Y.*P'); %%World GDP

%%%%%%%%%%%%%%%%%%%%%%%%%
%%With trade balance equation as the updating rule solve for wages (update wages)
%%%%%%%%%%%%%%%%%%%%%%%%%

%%First obtain royalties 

PiT= (1/(sig-1))*w.*LT; %%Profits T sector
PiNT= (1/(sig-1))*w.*LNT; %%Profits NT sector

Pi=(1/(sig-1))*w.*L; %%Total profits


r=(1+gg)^(1/(sig-1))/beta;  %%Discount factor (real interest rate)

for i=1:M
B(i)=0;
end



%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%
 %Endogenous chi


for i=1:M
     for j=1:M
         
chi(i,j) =rho(i);
     
     end
 end
 %%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%



%%Value of adopted technologies

for i=1:M
    for j=1:M
vadopter(i,j)=(1-chi(i,j))*(1-(1/r)*(1/(1))*((1+gg)^(1/(sig-1))/(1+gg)))^(-1)*PiT(i); %%Value adopted technology to adopter
vinnovator(i,j)=(chi(i,j))*(1-(1/r)*(1/(1))*((1+gg)^(1/(sig-1))/(1+gg)))^(-1)*PiT(i); %%Value adopted technology innovator
    end
end

 %%Guess H_{in}^a 

diffa=10; itera=0; 

     
%...GIT'R DUN!...
 %display('---------------------------------------------------------------');
 
while diffa>10^-6;
     itera=itera+1;
    
%%Probability of adoption and value of an unadopte dtechnology

for i=1:M
    for j=1:M
   epsbar(i,j) = eps(i,j)/((P(i)*Ha(i,j)/Yw)^betaa);     
   eps_test(i,j)=epsbar(i,j)*(P(i)*Ha(i,j)/Yw)^betaa;     %%Probability of adoption
   Jadopter(i,j) = (1-(1/r)*(1/(1))*((1+gg)^(1/(sig-1))/(1+gg))*(eps(i,j)*betaa+(1-eps(i,j))))^(-1)*(eps(i,j)*(1/r)*(1/(1))*((1+gg)^(1/(sig-1))/(1+gg))*(1-betaa))*vadopter(i,j); %%Value not yet adopted technology adopter
   Jinnovator(i,j) = (1-(1/r)*(1/(1))*((1+gg)^(1/(sig-1))/(1+gg))*((1-eps(i,j))))^(-1)*(eps(i,j)*(1/r)*(1/(1))*((1+gg)^(1/(sig-1))/(1+gg)))*vinnovator(i,j); %%Value not-yet adopted technology innovator (innovator chooses Hr based on this)
   
    end
end

 
%%Value of innovation

for i=1:M
    for j=1:M
   vinnov(i,j)=Jinnovator(i,j)*(Tnew(j)/Tnew(i)); %%Value of an innovation (value of a patent)
    end
end
    
vinnov=sum(vinnov)';

for i=1:M
H_r(i) = Hrdata(i)*Y(i);
lam(i) = (P(i)*H_r(i))^(1-betar)*(vinnov(i)*betar*(Yw)^(-betar))^(-1);
H_r_test(i)=P(i)^-1*(vinnov(i)*betar*lam(i)*(Yw)^(-betar))^(1/(1-betar)); %%FOC R&D spending
end


%%Get Z
ZZ=(lam'./gg).*(P'.*H_r'./(Yw)).^betar.*Tnew; %%Growth newly invente dvarieties

%%Get A_Z

for i=1:M
    for j=1:M
A_Z1(i,j)=((eps(i,j))*ZZ(j)*(1+gg))/(eps(i,j)+gg); %%A(i,j)
    end
end

%A_Z2=sum(A_Z1');

for i=1:M
    for j=1:M
A_Z(i,j)=A_Z1(i,j)/Tnew(i); %%A(i,j)/T(i)
    end
end





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Update H_in^a
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%Use FOC of adoption to update adoption spending

for i=1:M
    for j=1:M
       LHS_Ha(i,j)=Ha(i,j)*P(i);
    end
end
RHS_Ha=gg*betaa.*(A_Z).*(1/r).*(vadopter-Jadopter).*(1/(1)).*((1+gg)^(1/(sig-1))/(1+gg));



Za=LHS_Ha-RHS_Ha;

%% If excess < 0 country one is importing less, the increase wage in country 1 and decrease wage in 2

 % Update wages
    scl = 0.5; % Scale factor used to ensure that Tw > 0.
    
    THa = Ha .* (1 - scl.*Za); 
    
    % Compute distance between 'new' and 'old' wage vector.
   
     diffa = (1/scl).*max(max((abs((THa-Ha)./THa))));
   
    
    % Update wage using smoothing parameter (convex combo b/w new and old).
    smooth =  .8;
    Ha = smooth.*THa + (1-smooth).*Ha;
    
end


%%%Update wages using trade balance equation (exports T sector + exports of
%%%T sector + royalties recieves= Imports T sector + Imports NT sector +
%%%royalties paid)

  LHStemp(1:M,1) = 0;
  Lroytemp(1:M,1:M)=0;
  RHStemp(1:M,1) = 0;
  Rroytemp(1:M,1:M)=0;
  IM_T(1:M,1) = 0;
  EX_T(1:M,1) = 0;
  IM_NT(1:M,1) = 0;
  EX_NT(1:M,1) = 0;
  IBT(1:M,1) = 0;

for i=1:M
    for j=1:M
        A(i,j)=A_Z1(i,j);
        roy(i,j) =chi(i,j)*ipr(i)*PiT(i)*A_Z(i,j); %%Royalties (fraction of profits each period) Think taht innovator sells patent to someone that then licenses the product
    end
end

for i=1:M
            for j =1:M
                if j~=i
                        IM_T(i)= IM_T(i)+ pi_shareT(i,j)*P(i)*YT(i)/(1+tau(i,j));
                         IBT(i)=IBT(i)+tau(i,j)*pi_shareT(i,j)*P(i)*YT(i)/(1+tau(i,j));
                        LHStemp(i)= LHStemp(i) +  pi_shareT(i,j)*P(i)*YT(i)/(1+tau(i,j));
                        Lroytemp(i)= Lroytemp(i) + roy(i,j);
                        LHS(i)= LHStemp(i)+Lroytemp(i);
                end
            end
        end

 for i=1:M
            for j =1:M
                if j~=i
                        EX_T(i)= EX_T(i)+ pi_shareT(j,i)*P(j)*YT(j)/(1+tau(j,i));
                        RHStemp(i)= RHStemp(i) + pi_shareT(j,i)*P(j)*YT(j)/(1+tau(j,i));
                        Rroytemp(i)= Rroytemp(i) + roy(j,i);
                        RHS(i)= RHStemp(i)+Rroytemp(i);
                end
            end
        end

%% Excess demand


Z = [(LHS'-RHS')]';

%% If excess < 0 country one is importing less, the increase wage in country 1 and decrease wage in 2

 % Update wages
    scl = 5*10^(-2); % Scale factor used to ensure that Tw > 0.
    
    Tw = w .* (1 - scl*Z'); 
    
    % Compute distance between 'new' and 'old' wage vector.
    
    diff = 1/scl*max(abs((Tw-w)./w) );
    
    % Update wage using smoothing parameter (convex combo b/w new and old).
    smooth =  0.8;
    w = smooth.*Tw + (1-smooth).*w;
    w=w./w(M);
     
    if mod(iter,100)==0
        telapsed=toc; sec=mod(telapsed,60); min=floor(telapsed/60);
        display(sprintf('\t\t SS iterations on w completed: %6.0f',iter));
        display(sprintf('\t\t\t time elapsed: %6.0f min, %2.0f sec',min,sec));
        display(sprintf('\t\t\t dif: %3.10f',diff));
    end
    
end

 %%Consumption
    
   C=Y-H_r'./P'-sum(Ha')';